home *** CD-ROM | disk | FTP | other *** search
/ Amiga Plus 2000 #5 / Amiga Plus CD - 2000 - No. 5.iso / Tools / Dev / lame_src / psymodel.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-01-01  |  37.6 KB  |  1,442 lines
  1. /*
  2.  *    psymodel.c
  3.  *
  4.  *    Copyright (c) 1999 Mark Taylor
  5.  *
  6.  * This library is free software; you can redistribute it and/or
  7.  * modify it under the terms of the GNU Library General Public
  8.  * License as published by the Free Software Foundation; either
  9.  * version 2 of the License, or (at your option) any later version.
  10.  *
  11.  * This library is distributed in the hope that it will be useful,
  12.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.     See the GNU
  14.  * Library General Public License for more details.
  15.  *
  16.  * You should have received a copy of the GNU Library General Public
  17.  * License along with this library; if not, write to the
  18.  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  19.  * Boston, MA 02111-1307, USA.
  20.  */
  21.  
  22. #include "util.h"
  23. #include "encoder.h"
  24. #include "psymodel.h"
  25. #include "l3side.h"
  26. #include <assert.h>
  27. #include "tables.h"
  28. #include "fft.h"
  29.  
  30. #ifdef M_LN10
  31. #define        LN_TO_LOG10        (M_LN10/10)
  32. #else
  33. #define         LN_TO_LOG10             0.2302585093
  34. #endif
  35.  
  36.  
  37. int L3para_read( lame_global_flags *gfp,
  38.           FLOAT8 sfreq, int numlines[CBANDS],int numlines_s[CBANDS], 
  39.           FLOAT8 minval[CBANDS], 
  40.           FLOAT8 s3_l[CBANDS][CBANDS],
  41.           FLOAT8 s3_s[CBANDS][CBANDS],
  42.           FLOAT8 SNR_s[CBANDS],
  43.           int bu_l[SBPSY_l], int bo_l[SBPSY_l],
  44.           FLOAT8 w1_l[SBPSY_l], FLOAT8 w2_l[SBPSY_l],
  45.           int bu_s[SBPSY_s], int bo_s[SBPSY_s],
  46.           FLOAT8 w1_s[SBPSY_s], FLOAT8 w2_s[SBPSY_s],
  47.           int *, int *, int *, int *);
  48.                                     
  49.  
  50.  
  51. int L3psycho_anal( lame_global_flags *gfp,
  52.                     short int *buffer[2],int gr_out , 
  53.                     FLOAT8 *ms_ratio,
  54.                     FLOAT8 *ms_ratio_next,
  55.             FLOAT8 *ms_ener_ratio,
  56.             III_psy_ratio masking_ratio[2][2],
  57.             III_psy_ratio masking_MS_ratio[2][2],
  58.             FLOAT8 percep_entropy[2],FLOAT8 percep_MS_entropy[2], 
  59.                     int blocktype_d[2])
  60. {
  61.   lame_internal_flags *gfc=gfp->internal_flags;
  62.  
  63. /* to get a good cache performance, one has to think about
  64.  * the sequence, in which the variables are used.  
  65.  * (Note: these static variables have been moved to the gfc-> struct,
  66.  * and their order in memory is layed out in util.h)
  67.  */
  68.   
  69.  
  70.   /* fft and energy calculation   */
  71.   FLOAT (*wsamp_l)[BLKSIZE];
  72.   FLOAT (*wsamp_s)[3][BLKSIZE_s];
  73.   FLOAT tot_ener[4];
  74.  
  75.   /* convolution   */
  76.   FLOAT8 eb[CBANDS];
  77.   FLOAT8 cb[CBANDS];
  78.   FLOAT8 thr[CBANDS];
  79.   
  80.   /* ratios    */
  81.   FLOAT8 ms_ratio_l=0,ms_ratio_s=0;
  82.  
  83.   /* block type  */
  84.   int blocktype[2],uselongblock[2];
  85.   
  86.   /* usual variables like loop indices, etc..    */
  87.   int numchn, chn, samplerate;
  88.   int   b, i, j, k;
  89.   int    sb,sblock;
  90.   FLOAT cwlimit;
  91.  
  92.   /*  use a simplified spreading function: */
  93.     /*#define NEWS3  */ 
  94. #if 1
  95.     /* AAC values, results in more masking over MP3 values */
  96. # define TMN 18
  97. # define NMT 6
  98. #else
  99.     /* MP3 values */
  100. # define TMN 29
  101. # define NMT 6
  102. #endif
  103.  
  104.  
  105.  
  106.  
  107.  
  108.   if(gfc->psymodel_init==0) {
  109.     FLOAT8    SNR_s[CBANDS];
  110.     gfc->psymodel_init=1;
  111.  
  112.  
  113.     samplerate = gfp->out_samplerate;
  114.     switch(gfp->out_samplerate){
  115.     case 32000: break;
  116.     case 44100: break;
  117.     case 48000: break;
  118.     case 16000: break;
  119.     case 22050: break;
  120.     case 24000: break;
  121.     case  8000: samplerate *= 2; break;  /* kludge so mpeg2.5 uses mpeg2 tables  for now */
  122.     case 11025: samplerate *= 2; break;
  123.     case 12000: samplerate *= 2; break;
  124.     default:    ERRORF("error, invalid sampling frequency: %d Hz\n",
  125.             gfp->out_samplerate);
  126.     return -1;
  127.     }
  128.  
  129.     gfc->ms_ener_ratio_old=.25;
  130.     gfc->blocktype_old[0]=STOP_TYPE;
  131.     gfc->blocktype_old[1]=STOP_TYPE;
  132.     gfc->blocktype_old[0]=SHORT_TYPE;
  133.     gfc->blocktype_old[1]=SHORT_TYPE;
  134.  
  135.     for (i=0; i<4; ++i) {
  136.       for (j=0; j<CBANDS; ++j) {
  137.     gfc->nb_1[i][j]=1e20;
  138.     gfc->nb_2[i][j]=1e20;
  139.       }
  140.       for ( sb = 0; sb < SBPSY_l; sb++ ) {
  141.     gfc->en[i].l[sb] = 1e20;
  142.     gfc->thm[i].l[sb] = 1e20;
  143.       }
  144.       for (j=0; j<3; ++j) {
  145.     for ( sb = 0; sb < SBPSY_s; sb++ ) {
  146.       gfc->en[i].s[sb][j] = 1e20;
  147.       gfc->thm[i].s[sb][j] = 1e20;
  148.     }
  149.       }
  150.     }
  151.  
  152.  
  153.  
  154.     
  155.     /*  gfp->cwlimit = sfreq*j/1024.0;  */
  156.     gfc->cw_lower_index=6;
  157.     if (gfp->cwlimit>0) 
  158.       cwlimit=gfp->cwlimit;
  159.     else
  160.       cwlimit=8.8717;
  161.     gfc->cw_upper_index = cwlimit*1000.0*1024.0/((FLOAT8)samplerate);
  162.     gfc->cw_upper_index=Min(HBLKSIZE-4,gfc->cw_upper_index);      /* j+3 < HBLKSIZE-1 */
  163.     gfc->cw_upper_index=Max(6,gfc->cw_upper_index);
  164.  
  165.     for ( j = 0; j < HBLKSIZE; j++ )
  166.       gfc->cw[j] = 0.4;
  167.     
  168.     
  169.  
  170.     i=L3para_read( gfp,(FLOAT8) samplerate,gfc->numlines_l,gfc->numlines_s,
  171.           gfc->minval,gfc->s3_l,gfc->s3_s,SNR_s,gfc->bu_l,
  172.           gfc->bo_l,gfc->w1_l,gfc->w2_l, gfc->bu_s,gfc->bo_s,
  173.           gfc->w1_s,gfc->w2_s,&gfc->npart_l_orig,&gfc->npart_l,
  174.           &gfc->npart_s_orig,&gfc->npart_s );
  175.     if (i!=0) return -1;
  176.     
  177.  
  178.  
  179.     /* npart_l_orig   = number of partition bands before convolution */
  180.     /* npart_l  = number of partition bands after convolution */
  181.     
  182.     for (i=0; i<gfc->npart_l; i++) {
  183.       for (j = 0; j < gfc->npart_l_orig; j++) {
  184.     if (gfc->s3_l[i][j] != 0.0)
  185.       break;
  186.       }
  187.       gfc->s3ind[i][0] = j;
  188.       
  189.       for (j = gfc->npart_l_orig - 1; j > 0; j--) {
  190.     if (gfc->s3_l[i][j] != 0.0)
  191.       break;
  192.       }
  193.       gfc->s3ind[i][1] = j;
  194.     }
  195.  
  196.  
  197.     for (i=0; i<gfc->npart_s; i++) {
  198.       for (j = 0; j < gfc->npart_s_orig; j++) {
  199.     if (gfc->s3_s[i][j] != 0.0)
  200.       break;
  201.       }
  202.       gfc->s3ind_s[i][0] = j;
  203.       
  204.       for (j = gfc->npart_s_orig - 1; j > 0; j--) {
  205.     if (gfc->s3_s[i][j] != 0.0)
  206.       break;
  207.       }
  208.       gfc->s3ind_s[i][1] = j;
  209.     }
  210.     
  211.     
  212.     /*  
  213.       #include "debugscalefac.c"
  214.     */
  215.     
  216.  
  217.  
  218.  
  219. #define rpelev 2
  220. #define rpelev2 16
  221.  
  222.     /* compute norm_l, norm_s instead of relying on table data */
  223.     for ( b = 0;b < gfc->npart_l; b++ ) {
  224.       FLOAT8 norm=0;
  225.       for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ ) {
  226.     norm += gfc->s3_l[b][k];
  227.       }
  228.       for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ ) {
  229.     gfc->s3_l[b][k] /=  norm;
  230.       }
  231.       /*DEBUGF("%i  norm=%f  norm_l=%f \n",b,1/norm,norm_l[b]);*/
  232.     }
  233.  
  234.  
  235.     for ( b = 0;b < gfc->npart_s; b++ ) {
  236.       FLOAT8 norm=0;
  237.       for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ ) {
  238.     norm += gfc->s3_s[b][k];
  239.       }
  240.       for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ ) {
  241.     gfc->s3_s[b][k] *= SNR_s[b] / norm;
  242.       }
  243.       /*DEBUGF("%i  norm=%f  norm_s=%f \n",b,1/norm,norm_l[b]);*/
  244.     }
  245.     
  246.     init_fft();
  247.   }
  248.   /************************* End of Initialization *****************************/
  249.   
  250.  
  251.  
  252.   
  253.   
  254.   numchn = gfc->stereo;
  255.   /* chn=2 and 3 = Mid and Side channels */
  256.   if (gfp->mode == MPG_MD_JOINT_STEREO) numchn=4;
  257.   for (chn=0; chn<numchn; chn++) {
  258.  
  259.     /* there is a one granule delay.  Copy maskings computed last call
  260.      * into masking_ratio to return to calling program.
  261.      */
  262.     if (chn<2) {    
  263.       /* LR maskings  */
  264.       percep_entropy[chn] = gfc->pe[chn]; 
  265.       masking_ratio[gr_out][chn].thm = gfc->thm[chn];
  266.       masking_ratio[gr_out][chn].en = gfc->en[chn];
  267.     }else{
  268.       /* MS maskings  */
  269.       percep_MS_entropy[chn-2] = gfc->pe[chn]; 
  270.       masking_MS_ratio[gr_out][chn-2].en = gfc->en[chn];
  271.       masking_MS_ratio[gr_out][chn-2].thm = gfc->thm[chn];
  272.     }
  273.       
  274.  
  275.     /**********************************************************************
  276.      *  compute FFTs
  277.      **********************************************************************/
  278.     wsamp_s = gfc->wsamp_S+(chn & 1);
  279.     wsamp_l = gfc->wsamp_L+(chn & 1);
  280.     if (chn<2) {    
  281.       fft_long ( *wsamp_l, chn, buffer);
  282.       fft_short( *wsamp_s, chn, buffer);
  283.     } 
  284.     /* FFT data for mid and side channel is derived from L & R */
  285.     if (chn == 2)
  286.       {
  287.         for (j = BLKSIZE-1; j >=0 ; --j)
  288.         {
  289.           FLOAT l = gfc->wsamp_L[0][j];
  290.           FLOAT r = gfc->wsamp_L[1][j];
  291.           gfc->wsamp_L[0][j] = (l+r)*(FLOAT)(SQRT2*0.5);
  292.           gfc->wsamp_L[1][j] = (l-r)*(FLOAT)(SQRT2*0.5);
  293.         }
  294.         for (b = 2; b >= 0; --b)
  295.         {
  296.           for (j = BLKSIZE_s-1; j >= 0 ; --j)
  297.           {
  298.             FLOAT l = gfc->wsamp_S[0][b][j];
  299.             FLOAT r = gfc->wsamp_S[1][b][j];
  300.             gfc->wsamp_S[0][b][j] = (l+r)*(FLOAT)(SQRT2*0.5);
  301.             gfc->wsamp_S[1][b][j] = (l-r)*(FLOAT)(SQRT2*0.5);
  302.           }
  303.         }
  304.       }
  305.  
  306.  
  307.     /**********************************************************************
  308.      *  compute energies
  309.      **********************************************************************/
  310.     
  311.     
  312.     
  313.     gfc->energy[0]  = (*wsamp_l)[0];
  314.     gfc->energy[0] *= gfc->energy[0];
  315.     
  316.     tot_ener[chn] = gfc->energy[0]; /* sum total energy at nearly no extra cost */
  317.     
  318.     for (j=BLKSIZE/2-1; j >= 0; --j)
  319.     {
  320.       FLOAT re = (*wsamp_l)[BLKSIZE/2-j];
  321.       FLOAT im = (*wsamp_l)[BLKSIZE/2+j];
  322.       gfc->energy[BLKSIZE/2-j] = (re * re + im * im) * (FLOAT)0.5;
  323.  
  324.       if (BLKSIZE/2-j > 10)
  325.     tot_ener[chn] += gfc->energy[BLKSIZE/2-j];
  326.     }
  327.     for (b = 2; b >= 0; --b)
  328.     {
  329.       gfc->energy_s[b][0]  = (*wsamp_s)[b][0];
  330.       gfc->energy_s[b][0] *=  gfc->energy_s [b][0];
  331.       for (j=BLKSIZE_s/2-1; j >= 0; --j)
  332.       {
  333.         FLOAT re = (*wsamp_s)[b][BLKSIZE_s/2-j];
  334.         FLOAT im = (*wsamp_s)[b][BLKSIZE_s/2+j];
  335.         gfc->energy_s[b][BLKSIZE_s/2-j] = (re * re + im * im) * (FLOAT)0.5;
  336.       }
  337.     }
  338.  
  339.  
  340.   if(gfp->gtkflag) {
  341.     for (j=0; j<HBLKSIZE ; j++) {
  342.       gfc->pinfo->energy[gr_out][chn][j]=gfc->energy_save[chn][j];
  343.       gfc->energy_save[chn][j]=gfc->energy[j];
  344.     }
  345.   }
  346.     
  347.     /**********************************************************************
  348.      *    compute unpredicatability of first six spectral lines            * 
  349.      **********************************************************************/
  350.     for ( j = 0; j < gfc->cw_lower_index; j++ )
  351.       {     /* calculate unpredictability measure cw */
  352.     FLOAT an, a1, a2;
  353.     FLOAT bn, b1, b2;
  354.     FLOAT rn, r1, r2;
  355.     FLOAT numre, numim, den;
  356.  
  357.     a2 = gfc-> ax_sav[chn][1][j];
  358.     b2 = gfc-> bx_sav[chn][1][j];
  359.     r2 = gfc-> rx_sav[chn][1][j];
  360.     a1 = gfc-> ax_sav[chn][1][j] = gfc-> ax_sav[chn][0][j];
  361.     b1 = gfc-> bx_sav[chn][1][j] = gfc-> bx_sav[chn][0][j];
  362.     r1 = gfc-> rx_sav[chn][1][j] = gfc-> rx_sav[chn][0][j];
  363.     an = gfc-> ax_sav[chn][0][j] = (*wsamp_l)[j];
  364.     bn = gfc-> bx_sav[chn][0][j] = j==0 ? (*wsamp_l)[0] : (*wsamp_l)[BLKSIZE-j];  
  365.     rn = gfc-> rx_sav[chn][0][j] = sqrt(gfc->energy[j]);
  366.  
  367.     { /* square (x1,y1) */
  368.       if( r1 != 0 ) {
  369.         numre = (a1*b1);
  370.         numim = (a1*a1-b1*b1)*(FLOAT)0.5;
  371.         den = r1*r1;
  372.       } else {
  373.         numre = 1;
  374.         numim = 0;
  375.         den = 1;
  376.       }
  377.     }
  378.     
  379.     { /* multiply by (x2,-y2) */
  380.       if( r2 != 0 ) {
  381.         FLOAT tmp2 = (numim+numre)*(a2+b2)*(FLOAT)0.5;
  382.         FLOAT tmp1 = -a2*numre+tmp2;
  383.         numre =       -b2*numim+tmp2;
  384.         numim = tmp1;
  385.         den *= r2;
  386.       } else {
  387.         /* do nothing */
  388.       }
  389.     }
  390.     
  391.     { /* r-prime factor */
  392.       FLOAT tmp = (2*r1-r2)/den;
  393.       numre *= tmp;
  394.       numim *= tmp;
  395.     }
  396.     den=rn+fabs(2*r1-r2);
  397.     if( den != 0 ) {
  398.       numre = (an+bn)*(FLOAT)0.5-numre;
  399.       numim = (an-bn)*(FLOAT)0.5-numim;
  400.       den = sqrt(numre*numre+numim*numim)/den;
  401.     }
  402.     gfc->cw[j] = den;
  403.       }
  404.  
  405.  
  406.  
  407.     /**********************************************************************
  408.      *     compute unpredicatibility of next 200 spectral lines            *
  409.      **********************************************************************/ 
  410.     for ( j = gfc->cw_lower_index; j < gfc->cw_upper_index; j += 4 )
  411.       {/* calculate unpredictability measure cw */
  412.     FLOAT rn, r1, r2;
  413.     FLOAT numre, numim, den;
  414.     
  415.     k = (j+2) / 4; 
  416.     
  417.     { /* square (x1,y1) */
  418.       r1 = gfc->energy_s[0][k];
  419.       if( r1 != 0 ) {
  420.         FLOAT a1 = (*wsamp_s)[0][k]; 
  421.         FLOAT b1 = (*wsamp_s)[0][BLKSIZE_s-k]; /* k is never 0 */
  422.         numre = (a1*b1);
  423.         numim = (a1*a1-b1*b1)*(FLOAT)0.5;
  424.         den = r1;
  425.         r1 = sqrt(r1);
  426.       } else {
  427.         numre = 1;
  428.         numim = 0;
  429.         den = 1;
  430.       }
  431.     }
  432.     
  433.     
  434.     { /* multiply by (x2,-y2) */
  435.       r2 = gfc->energy_s[2][k];
  436.       if( r2 != 0 ) {
  437.         FLOAT a2 = (*wsamp_s)[2][k]; 
  438.         FLOAT b2 = (*wsamp_s)[2][BLKSIZE_s-k];
  439.         
  440.         
  441.         FLOAT tmp2 = (numim+numre)*(a2+b2)*(FLOAT)0.5;
  442.         FLOAT tmp1 = -a2*numre+tmp2;
  443.         numre =       -b2*numim+tmp2;
  444.         numim = tmp1;
  445.         
  446.         r2 = sqrt(r2);
  447.         den *= r2;
  448.       } else {
  449.         /* do nothing */
  450.       }
  451.     }
  452.     
  453.     { /* r-prime factor */
  454.       FLOAT tmp = (2*r1-r2)/den;
  455.       numre *= tmp;
  456.       numim *= tmp;
  457.     }
  458.     
  459.     rn = sqrt(gfc->energy_s[1][k]);
  460.     den=rn+fabs(2*r1-r2);
  461.     if( den != 0 ) {
  462.       FLOAT an = (*wsamp_s)[1][k]; 
  463.       FLOAT bn = (*wsamp_s)[1][BLKSIZE_s-k];
  464.       numre = (an+bn)*(FLOAT)0.5-numre;
  465.       numim = (an-bn)*(FLOAT)0.5-numim;
  466.       den = sqrt(numre*numre+numim*numim)/den;
  467.     }
  468.     
  469.     gfc->cw[j+1] = gfc->cw[j+2] = gfc->cw[j+3] = gfc->cw[j] = den;
  470.       }
  471.     
  472. #if 0
  473.     for ( j = 14; j < HBLKSIZE-4; j += 4 )
  474.       {/* calculate energy from short ffts */
  475.     FLOAT8 tot,ave;
  476.     k = (j+2) / 4; 
  477.     for (tot=0, sblock=0; sblock < 3; sblock++)
  478.       tot+=gfc->energy_s[sblock][k];
  479.     ave = gfc->energy[j+1]+ gfc->energy[j+2]+ gfc->energy[j+3]+ gfc->energy[j];
  480.     ave /= 4.;
  481.     /*
  482.       DEBUGF("energy / tot %i %5.2f   %e  %e\n",j,ave/(tot*16./3.),
  483.       ave,tot*16./3.);
  484.     */
  485.     gfc->energy[j+1] = gfc->energy[j+2] = gfc->energy[j+3] =  gfc->energy[j]=tot;
  486.       }
  487. #endif
  488.     
  489.     
  490.     
  491.     
  492.     
  493.     
  494.     
  495.     
  496.     /**********************************************************************
  497.      *    Calculate the energy and the unpredictability in the threshold   *
  498.      *    calculation partitions                                           *
  499.      **********************************************************************/
  500.     b = 0;
  501.     for (j = 0; j < gfc->cw_upper_index && gfc->numlines_l[b] && b<gfc->npart_l_orig; )
  502.       {
  503.     FLOAT8 ebb, cbb;
  504.  
  505.     ebb = gfc->energy[j];
  506.     cbb = gfc->energy[j] * gfc->cw[j];
  507.     j++;
  508.  
  509.  
  510.     for (i = gfc->numlines_l[b] - 1; i > 0; i--)
  511.       {
  512.         ebb += gfc->energy[j];
  513.         cbb += gfc->energy[j] * gfc->cw[j];
  514.         j++;
  515.       }
  516.     eb[b] = ebb;
  517.     cb[b] = cbb;
  518.     b++;
  519.       }
  520.  
  521.     for (; b < gfc->npart_l_orig; b++ )
  522.       {
  523.     FLOAT8 ebb = gfc->energy[j++];
  524.     assert(gfc->numlines_l[b]);
  525.     for (i = gfc->numlines_l[b] - 1; i > 0; i--)
  526.       {
  527.         ebb += gfc->energy[j++];
  528.       }
  529.     eb[b] = ebb;
  530.     cb[b] = ebb * 0.4;
  531.       }
  532.  
  533.  
  534.     /**********************************************************************
  535.      *      convolve the partitioned energy and unpredictability           *
  536.      *      with the spreading function, s3_l[b][k]                        *
  537.      ******************************************************************** */
  538.     gfc->pe[chn] = 0;        /*  calculate percetual entropy */
  539.     for ( b = 0;b < gfc->npart_l; b++ )
  540.       {
  541.     FLOAT8 tbb,ecb,ctb;
  542.  
  543.     ecb = 0;
  544.     ctb = 0;
  545.     for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ )
  546.       {
  547.         ecb += gfc->s3_l[b][k] * eb[k];    /* sprdngf for Layer III */
  548.         ctb += gfc->s3_l[b][k] * cb[k];
  549.       }
  550.  
  551.     /* calculate the tonality of each threshold calculation partition 
  552.      * calculate the SNR in each threshhold calculation partition 
  553.      * tonality = -0.299 - .43*log(ctb/ecb);
  554.      * tonality = 0:           use NMT   (lots of masking)
  555.      * tonality = 1:           use TMN   (little masking)
  556.      */
  557.  
  558. /* ISO values */
  559. #define CONV1 (-.299)
  560. #define CONV2 (-.43)
  561.  
  562.     tbb = ecb;
  563.     if (tbb != 0)
  564.       {
  565.         tbb = ctb / tbb;
  566.         if (tbb <= exp((1-CONV1)/CONV2)) 
  567.           { /* tonality near 1 */
  568.         tbb = exp(-LN_TO_LOG10 * TMN);
  569.           }
  570.         else if (tbb >= exp((0-CONV1)/CONV2)) 
  571.           {/* tonality near 0 */
  572.         tbb = exp(-LN_TO_LOG10 * NMT);
  573.           }
  574.         else
  575.           {
  576.         /* convert to tonality index */
  577.         /* tonality small:   tbb=1 */
  578.         /* tonality large:   tbb=-.299 */
  579.         tbb = CONV1 + CONV2*log(tbb);
  580.         tbb = exp(-LN_TO_LOG10 * ( TMN*tbb + (1-tbb)*NMT) );
  581.           }
  582.       }
  583.     /* at this point, tbb represents the amount the spreading function
  584.      * will be reduced.  The smaller the value, the less masking.
  585.      * minval[] = 1 (0db)     says just use tbb.
  586.      * minval[]= .01 (-20db)  says reduce spreading function by 
  587.      *                        at least 20db.  
  588.      */
  589.     tbb = Min(gfc->minval[b], tbb);
  590.     ecb *= tbb;
  591.  
  592.  
  593.  
  594.     /* long block pre-echo control.   */
  595.     /* rpelev=2.0, rpelev2=16.0 */
  596.     /* note: all surges in PE are because of this pre-echo formula
  597.      * for thr[b].  If it this is not used, PE is always around 600
  598.      */
  599.     /* dont use long block pre-echo control if previous granule was 
  600.      * a short block.  This is to avoid the situation:   
  601.      * frame0:  quiet (very low masking)  
  602.      * frame1:  surge  (triggers short blocks)
  603.      * frame2:  regular frame.  looks like pre-echo when compared to 
  604.      *          frame0, but all pre-echo was in frame1.
  605.      */
  606.     //    ecb = Max(ecb,gfc->ATH_partitionbands[b]);
  607.     if (gfc->blocktype_old[chn] == SHORT_TYPE )
  608.       thr[b] = ecb; /* Min(ecb, rpelev*gfc->nb_1[chn][b]); */
  609.     else
  610.       thr[b] = Min(ecb, Min(rpelev*gfc->nb_1[chn][b],rpelev2*gfc->nb_2[chn][b]) );
  611.  
  612.     gfc->nb_2[chn][b] = gfc->nb_1[chn][b];
  613.     gfc->nb_1[chn][b] = ecb;
  614.  
  615.     {
  616.       FLOAT8 thrpe;
  617.       thrpe = Max(thr[b],gfc->ATH_partitionbands[b]);
  618.       /*
  619.         printf("%i thr=%e   ATH=%e  \n",b,thr[b],gfc->ATH_partitionbands[b]);
  620.       */
  621.       if (thrpe < eb[b])
  622.         gfc->pe[chn] -= gfc->numlines_l[b] * log(thrpe / eb[b]);
  623.     }
  624.       }
  625.  
  626.  
  627.     
  628.     /*************************************************************** 
  629.      * determine the block type (window type) based on L & R channels
  630.      * 
  631.      ***************************************************************/
  632.     {  /* compute PE for all 4 channels */
  633.       FLOAT mn,mx,ma=0,mb=0,mc=0;
  634.       
  635.       for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++)
  636.     {
  637.       ma += gfc->energy_s[0][j];
  638.       mb += gfc->energy_s[1][j];
  639.       mc += gfc->energy_s[2][j];
  640.     }
  641.       mn = Min(ma,mb);
  642.       mn = Min(mn,mc);
  643.       mx = Max(ma,mb);
  644.       mx = Max(mx,mc);
  645.       
  646.       /* bit allocation is based on pe.  */
  647.       if (mx>mn) {
  648.     FLOAT8 tmp = 400*log(mx/(1e-12+mn));
  649.     if (tmp>gfc->pe[chn]) gfc->pe[chn]=tmp;
  650.       }
  651.  
  652.       /* block type is based just on L or R channel */      
  653.       if (chn<2) {
  654.     uselongblock[chn] = 1;
  655.     
  656.     /* tuned for t1.wav.  doesnt effect most other samples */
  657.     if (gfc->pe[chn] > 3000) 
  658.       uselongblock[chn]=0;
  659.     
  660.     if ( mx > 30*mn ) 
  661.       {/* big surge of energy - always use short blocks */
  662.         uselongblock[chn] = 0;
  663.       } 
  664.     else if ((mx > 10*mn) && (gfc->pe[chn] > 1000))
  665.       {/* medium surge, medium pe - use short blocks */
  666.         uselongblock[chn] = 0;
  667.       }
  668.     
  669.     /* disable short blocks */
  670.     if (gfp->no_short_blocks)
  671.       uselongblock[chn]=1;
  672.       }
  673.     }
  674.  
  675.  
  676.     if (gfp->gtkflag) {
  677.       FLOAT mn,mx,ma=0,mb=0,mc=0;
  678.  
  679.       for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++)
  680.       {
  681.         ma += gfc->energy_s[0][j];
  682.         mb += gfc->energy_s[1][j];
  683.         mc += gfc->energy_s[2][j];
  684.       }
  685.       mn = Min(ma,mb);
  686.       mn = Min(mn,mc);
  687.       mx = Max(ma,mb);
  688.       mx = Max(mx,mc);
  689.  
  690.       gfc->pinfo->ers[gr_out][chn]=gfc->ers_save[chn];
  691.       gfc->ers_save[chn]=(mx/(1e-12+mn));
  692.       gfc->pinfo->pe[gr_out][chn]=gfc->pe_save[chn];
  693.       gfc->pe_save[chn]=gfc->pe[chn];
  694.     }
  695.  
  696.  
  697.     /*************************************************************** 
  698.      * compute masking thresholds for both short and long blocks
  699.      ***************************************************************/
  700.     /* longblock threshold calculation (part 2) */
  701.     for ( sb = 0; sb < SBPSY_l; sb++ )
  702.       {
  703. #if 1  /* additive masking */
  704.     FLOAT8 enn = gfc->w1_l[sb] * eb[gfc->bu_l[sb]] + gfc->w2_l[sb] * eb[gfc->bo_l[sb]];
  705.     FLOAT8 thmm = gfc->w1_l[sb] *thr[gfc->bu_l[sb]] + gfc->w2_l[sb] * thr[gfc->bo_l[sb]];
  706.     for ( b = gfc->bu_l[sb]+1; b < gfc->bo_l[sb]; b++ )
  707.       {
  708.         enn  += eb[b];
  709.         thmm += thr[b];
  710.       }
  711. #else
  712.     FLOAT8 enn = gfc->w1_l[sb] * eb[gfc->bu_l[sb]] + gfc->w2_l[sb] * eb[gfc->bo_l[sb]];
  713.     FLOAT8 thmm = Min(thr[gfc->bu_l[sb]],thr[gfc->bo_l[sb]]);
  714.     for ( b = gfc->bu_l[sb]+1; b < gfc->bo_l[sb]; b++ )
  715.       {
  716.         enn  += eb[b];
  717.         thmm = Min(thr[b],thmm);
  718.       }
  719.     thmm*=(1+gfc->bo_l[sb]-gfc->bu_l[sb]);
  720. #endif
  721.     gfc->en[chn].l[sb] = enn;
  722.     gfc->thm[chn].l[sb] = thmm;
  723.       }
  724.     
  725.  
  726.     /* threshold calculation for short blocks */
  727.     for ( sblock = 0; sblock < 3; sblock++ )
  728.       {
  729.     j = 0;
  730.     for ( b = 0; b < gfc->npart_s_orig; b++ )
  731.       {
  732.         FLOAT ecb = gfc->energy_s[sblock][j++];
  733.         for (i = 1 ; i<gfc->numlines_s[b]; ++i)
  734.           {
  735.         ecb += gfc->energy_s[sblock][j++];
  736.           }
  737.         eb[b] = ecb;
  738.       }
  739.  
  740.     for ( b = 0; b < gfc->npart_s; b++ )
  741.       {
  742.         FLOAT8 ecb = 0;
  743.         for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ )
  744.           {
  745.         ecb += gfc->s3_s[b][k] * eb[k];
  746.           }
  747.         thr[b] = Max (1e-6, ecb);
  748.       }
  749.  
  750.     for ( sb = 0; sb < SBPSY_s; sb++ )
  751.       {
  752.         FLOAT8 enn  = gfc->w1_s[sb] * eb[gfc->bu_s[sb]] + gfc->w2_s[sb] * eb[gfc->bo_s[sb]];
  753.         FLOAT8 thmm = gfc->w1_s[sb] *thr[gfc->bu_s[sb]] + gfc->w2_s[sb] * thr[gfc->bo_s[sb]];
  754.         for ( b = gfc->bu_s[sb]+1; b < gfc->bo_s[sb]; b++ )
  755.           {
  756.         enn  += eb[b];
  757.         thmm += thr[b];
  758.           }
  759.         gfc->en[chn].s[sb][sblock] = enn;
  760.         gfc->thm[chn].s[sb][sblock] = thmm;
  761.       }
  762.       }
  763.  
  764.   } /* end loop over chn */
  765.  
  766.  
  767.   /* compute M/S thresholds from Johnston & Ferreira 1992 ICASSP paper */
  768.   if ( numchn==4 /* mid/side and r/l */) {
  769.     FLOAT8 rside,rmid,mld;
  770.     int chmid=2,chside=3; 
  771.     
  772.     for ( sb = 0; sb < SBPSY_l; sb++ ) {
  773.       /* use this fix if L & R masking differs by 2db or less */
  774.       /* if db = 10*log10(x2/x1) < 2 */
  775.       /* if (x2 < 1.58*x1) { */
  776.       if (gfc->thm[0].l[sb] <= 1.58*gfc->thm[1].l[sb]
  777.       && gfc->thm[1].l[sb] <= 1.58*gfc->thm[0].l[sb]) {
  778.  
  779.     mld = gfc->mld_l[sb]*gfc->en[chside].l[sb];
  780.     rmid = Max(gfc->thm[chmid].l[sb], Min(gfc->thm[chside].l[sb],mld));
  781.  
  782.     mld = gfc->mld_l[sb]*gfc->en[chmid].l[sb];
  783.     rside = Max(gfc->thm[chside].l[sb],Min(gfc->thm[chmid].l[sb],mld));
  784.  
  785.     gfc->thm[chmid].l[sb]=rmid;
  786.     gfc->thm[chside].l[sb]=rside;
  787.       }
  788.     }
  789.     for ( sb = 0; sb < SBPSY_s; sb++ ) {
  790.       for ( sblock = 0; sblock < 3; sblock++ ) {
  791.     if (gfc->thm[0].s[sb][sblock] <= 1.58*gfc->thm[1].s[sb][sblock]
  792.         && gfc->thm[1].s[sb][sblock] <= 1.58*gfc->thm[0].s[sb][sblock]) {
  793.  
  794.       mld = gfc->mld_s[sb]*gfc->en[chside].s[sb][sblock];
  795.       rmid = Max(gfc->thm[chmid].s[sb][sblock],Min(gfc->thm[chside].s[sb][sblock],mld));
  796.  
  797.       mld = gfc->mld_s[sb]*gfc->en[chmid].s[sb][sblock];
  798.       rside = Max(gfc->thm[chside].s[sb][sblock],Min(gfc->thm[chmid].s[sb][sblock],mld));
  799.  
  800.       gfc->thm[chmid].s[sb][sblock]=rmid;
  801.       gfc->thm[chside].s[sb][sblock]=rside;
  802.     }
  803.       }
  804.     }
  805.   }
  806.  
  807.  
  808.   
  809.  
  810.   
  811.   
  812.   if (gfp->mode == MPG_MD_JOINT_STEREO)  {
  813.     /* determin ms_ratio from masking thresholds*/
  814.     /* use ms_stereo (ms_ratio < .35) if average thresh. diff < 5 db */
  815.     FLOAT8 db,x1,x2,sidetot=0,tot=0;
  816.     for (sb= SBPSY_l/4 ; sb< SBPSY_l; sb ++ ) {
  817.       x1 = Min(gfc->thm[0].l[sb],gfc->thm[1].l[sb]);
  818.       x2 = Max(gfc->thm[0].l[sb],gfc->thm[1].l[sb]);
  819.       /* thresholds difference in db */
  820.       if (x2 >= 1000*x1)  db=3;
  821.       else db = log10(x2/x1);  
  822.       /*  DEBUGF("db = %f %e %e  \n",db,gfc->thm[0].l[sb],gfc->thm[1].l[sb]);*/
  823.       sidetot += db;
  824.       tot++;
  825.     }
  826.     ms_ratio_l= (sidetot/tot)*0.7; /* was .35*(sidetot/tot)/5.0*10 */
  827.     ms_ratio_l = Min(ms_ratio_l,0.5);
  828.     
  829.     sidetot=0; tot=0;
  830.     for ( sblock = 0; sblock < 3; sblock++ )
  831.       for ( sb = SBPSY_s/4; sb < SBPSY_s; sb++ ) {
  832.     x1 = Min(gfc->thm[0].s[sb][sblock],gfc->thm[1].s[sb][sblock]);
  833.     x2 = Max(gfc->thm[0].s[sb][sblock],gfc->thm[1].s[sb][sblock]);
  834.     /* thresholds difference in db */
  835.     if (x2 >= 1000*x1)  db=3;
  836.     else db = log10(x2/x1);  
  837.     sidetot += db;
  838.     tot++;
  839.       }
  840.     ms_ratio_s = (sidetot/tot)*0.7; /* was .35*(sidetot/tot)/5.0*10 */
  841.     ms_ratio_s = Min(ms_ratio_s,.5);
  842.   }
  843.  
  844.   /*************************************************************** 
  845.    * determin final block type
  846.    ***************************************************************/
  847.  
  848.   for (chn=0; chn<gfc->stereo; chn++) {
  849.     blocktype[chn] = NORM_TYPE;
  850.   }
  851.  
  852.  
  853.   if (gfc->stereo==2) {
  854.     if (!gfp->allow_diff_short || gfp->mode==MPG_MD_JOINT_STEREO) {
  855.       /* force both channels to use the same block type */
  856.       /* this is necessary if the frame is to be encoded in ms_stereo.  */
  857.       /* But even without ms_stereo, FhG  does this */
  858.       int bothlong= (uselongblock[0] && uselongblock[1]);
  859.       if (!bothlong) {
  860.     uselongblock[0]=0;
  861.     uselongblock[1]=0;
  862.       }
  863.     }
  864.   }
  865.  
  866.   
  867.   
  868.   /* update the blocktype of the previous granule, since it depends on what
  869.    * happend in this granule */
  870.   for (chn=0; chn<gfc->stereo; chn++) {
  871.     if ( uselongblock[chn])
  872.       {                /* no attack : use long blocks */
  873.     switch( gfc->blocktype_old[chn] ) 
  874.       {
  875.       case NORM_TYPE:
  876.       case STOP_TYPE:
  877.         blocktype[chn] = NORM_TYPE;
  878.         break;
  879.       case SHORT_TYPE:
  880.         blocktype[chn] = STOP_TYPE; 
  881.         break;
  882.       case START_TYPE:
  883.         ERRORF( "Error in block selecting\n" );
  884.         LAME_ERROR_EXIT();
  885.         break; /* problem */
  886.       }
  887.       } else   {
  888.     /* attack : use short blocks */
  889.     blocktype[chn] = SHORT_TYPE;
  890.     if ( gfc->blocktype_old[chn] == NORM_TYPE ) {
  891.       gfc->blocktype_old[chn] = START_TYPE;
  892.     }
  893.     if ( gfc->blocktype_old[chn] == STOP_TYPE ) {
  894.       gfc->blocktype_old[chn] = SHORT_TYPE ;
  895.     }
  896.       }
  897.     
  898.     blocktype_d[chn] = gfc->blocktype_old[chn];  /* value returned to calling program */
  899.     gfc->blocktype_old[chn] = blocktype[chn];    /* save for next call to l3psy_anal */
  900.   }
  901.   
  902.   if (blocktype_d[0]==2) 
  903.     *ms_ratio = gfc->ms_ratio_s_old;
  904.   else
  905.     *ms_ratio = gfc->ms_ratio_l_old;
  906.  
  907.   gfc->ms_ratio_s_old = ms_ratio_s;
  908.   gfc->ms_ratio_l_old = ms_ratio_l;
  909.  
  910.   /* we dont know the block type of this frame yet - assume long */
  911.   *ms_ratio_next = ms_ratio_l;
  912.  
  913.  
  914.  
  915.   /*********************************************************************/
  916.   /* compute side_energy / (side+mid)_energy */
  917.   /* 0 = no energy in side channel */
  918.   /* .5 = half of total energy in side channel */
  919.   /*********************************************************************/
  920.   if (numchn==4)  {
  921.     FLOAT tmp = tot_ener[3]+tot_ener[2];
  922.     *ms_ener_ratio = gfc->ms_ener_ratio_old;
  923.     gfc->ms_ener_ratio_old=0;
  924.     if (tmp>0) gfc->ms_ener_ratio_old=tot_ener[3]/tmp;
  925.   } else
  926.     /* we didn't compute ms_ener_ratios */
  927.     *ms_ener_ratio = 0;
  928.  
  929.   return 0;
  930. }
  931.  
  932.  
  933.  
  934.  
  935.  
  936. int L3para_read(lame_global_flags *gfp, FLOAT8 sfreq, int *numlines_l,int *numlines_s, 
  937. FLOAT8 *minval,
  938. FLOAT8 s3_l[CBANDS][CBANDS], FLOAT8 s3_s[CBANDS][CBANDS],
  939. FLOAT8 *SNR, 
  940. int *bu_l, int *bo_l, FLOAT8 *w1_l, FLOAT8 *w2_l, 
  941. int *bu_s, int *bo_s, FLOAT8 *w1_s, FLOAT8 *w2_s,
  942. int *npart_l_orig,int *npart_l,int *npart_s_orig,int *npart_s)
  943. {
  944.   lame_internal_flags *gfc=gfp->internal_flags;
  945.   FLOAT8 freq_tp;
  946.   FLOAT8 bval_l[CBANDS], bval_s[CBANDS];
  947.   int   cbmax=0, cbmax_tp;
  948.   FLOAT8 *p = psy_data;
  949.   int  sbmax ;
  950.   int  i,j,k2,loop;
  951.   int freq_scale=1;
  952.   int partition[HBLKSIZE];
  953.  
  954.   /******************************************************************/
  955.   /* Read long block data */
  956.   /******************************************************************/
  957.   for(loop=0;loop<6;loop++)
  958.     {
  959.       freq_tp = *p++;
  960.       cbmax_tp = (int) *p++;
  961.       cbmax_tp++;
  962.  
  963.       if (sfreq == freq_tp/freq_scale )
  964.     {
  965.       cbmax = cbmax_tp;
  966.       for(i=0,k2=0;i<cbmax_tp;i++)
  967.         {
  968.           j = (int) *p++;
  969.           numlines_l[i] = (int) *p++;
  970.           minval[i] = exp(-((*p++) ) * LN_TO_LOG10);
  971.           /* qthr_l[i] = *p++ */ p++;
  972.           /* norm_l[i] = *p++*/ p++;
  973.           /* bval_l[i] = *p++; */ p++;
  974.           if (j!=i)
  975.         {
  976.           ERRORF("1. please check \"psy_data\"");
  977.           return -1;
  978.         }
  979.         }
  980.     }
  981.       else
  982.     p += cbmax_tp * 6;
  983.     }
  984.  
  985.   *npart_l_orig = cbmax;
  986.  
  987.   /* Read short block data */
  988.   for(loop=0;loop<6;loop++)
  989.     {
  990.       freq_tp = *p++;
  991.       cbmax_tp = (int) *p++;
  992.       cbmax_tp++;
  993.       if (sfreq == freq_tp/freq_scale )
  994.     {
  995.       cbmax = cbmax_tp;
  996.       for(i=0,k2=0;i<cbmax_tp;i++)
  997.         {
  998.           j = (int) *p++;
  999.           numlines_s[i] = (int) *p++;
  1000.           /* qthr_s[i] = *p++*/  p++;         
  1001.           /* norm_s[i] =*p++ */ p++;         
  1002.           SNR[i] = *p++;            
  1003.           /* bval_s[i] = *p++ */ p++;
  1004.           if (j!=i)
  1005.         {
  1006.           ERRORF("3. please check \"psy_data\"");
  1007.           return -1;
  1008.         }
  1009.         }
  1010.     }
  1011.       else
  1012.     p += cbmax_tp * 6;
  1013.     }
  1014.   *npart_s_orig = cbmax;
  1015.  
  1016.   /* MPEG1 SNR_s data is given in db, convert to energy */
  1017.   if (gfp->version == 1) {
  1018.     for ( i = 0;i < *npart_s_orig; i++ ) {
  1019.       SNR[i]=exp( (FLOAT8) SNR[i] * LN_TO_LOG10 );
  1020.     }
  1021.   }
  1022.  
  1023.  
  1024.  
  1025.  
  1026.   /* Read long block data for converting threshold calculation 
  1027.      partitions to scale factor bands */
  1028.  
  1029.   for(loop=0;loop<6;loop++)
  1030.     {
  1031.       freq_tp = *p++;
  1032.       sbmax =  (int) *p++;
  1033.       sbmax++;
  1034.  
  1035.       if (sfreq == freq_tp/freq_scale)
  1036.     {
  1037.       for(i=0;i<sbmax;i++)
  1038.         {
  1039.           j = (int) *p++;
  1040.           p++;             
  1041.           bu_l[i] = (int) *p++;
  1042.           bo_l[i] = (int) *p++;
  1043.           w1_l[i] = (FLOAT8) *p++;
  1044.           w2_l[i] = (FLOAT8) *p++;
  1045.           if (j!=i)
  1046.         { ERRORF("30:please check \"psy_data\"\n");
  1047.         return -1;
  1048.         }
  1049.  
  1050.           if (i!=0)
  1051.         if ( (fabs(1.0-w1_l[i]-w2_l[i-1]) > 0.01 ) )
  1052.           {
  1053.             ERRORF("31l: please check \"psy_data.\"\n"
  1054.                            "w1,w2: %f %f \n",w1_l[i],w2_l[i-1]);
  1055.             return -1;
  1056.           }
  1057.         }
  1058.     }
  1059.       else
  1060.     p += sbmax * 6;
  1061.     }
  1062.  
  1063.   /* Read short block data for converting threshold calculation 
  1064.      partitions to scale factor bands */
  1065.  
  1066.   for(loop=0;loop<6;loop++)
  1067.     {
  1068.       freq_tp = *p++;
  1069.       sbmax = (int) *p++;
  1070.       sbmax++;
  1071.  
  1072.       if (sfreq == freq_tp/freq_scale)
  1073.     {
  1074.       for(i=0;i<sbmax;i++)
  1075.         {
  1076.           j = (int) *p++;
  1077.           p++;
  1078.           bu_s[i] = (int) *p++;
  1079.           bo_s[i] = (int) *p++;
  1080.           w1_s[i] = *p++;
  1081.           w2_s[i] = *p++;
  1082.           if (j!=i)
  1083.         { ERRORF("30:please check \"psy_data\"\n");
  1084.         return -1;
  1085.         }
  1086.  
  1087.           if (i!=0)
  1088.         if ( (fabs(1.0-w1_s[i]-w2_s[i-1]) > 0.01 ) )
  1089.           { 
  1090.                   ERRORF("31s: please check \"psy_data.\"\n"
  1091.                          "w1,w2: %f %f \n",w1_s[i],w2_s[i-1]);
  1092.           return -1;
  1093.           }
  1094.         }
  1095.     }
  1096.       else
  1097.     p += sbmax * 6;
  1098.     }
  1099.  
  1100.   /******************************************************************/
  1101.   /* done reading table data */
  1102.   /******************************************************************/
  1103.  
  1104.  
  1105.  
  1106.  
  1107. #undef NOTABLES
  1108. #ifdef NOTABLES
  1109.   /* compute numlines */
  1110.   j=0;
  1111.   for(i=0;i<CBANDS;i++)
  1112.     {
  1113.       FLOAT8 ji, bark1,bark2,delbark=.34;
  1114.       int k,j2;
  1115.  
  1116.       j2 = j;
  1117.       j2 = Min(j2,BLKSIZE/2);
  1118.       
  1119.       do {
  1120.     /* find smallest j2 >= j so that  (bark - bark_l[i-1]) < delbark */
  1121.     ji = j;
  1122.     bark1 = freq2bark(sfreq*ji/BLKSIZE);
  1123.     
  1124.     ++j2;
  1125.     ji = j2;
  1126.     bark2  = freq2bark(sfreq*ji/BLKSIZE);
  1127.       } while ((bark2 - bark1) < delbark  && j2<=BLKSIZE/2);
  1128.  
  1129.       /*
  1130.       DEBUGF("%i old n=%i  %f old numlines:  %i   new=%i (%i,%i) (%f,%f) \n",
  1131. i,*npart_l_orig,freq,numlines_l[i],j2-j,j,j2-1,bark1,bark2);
  1132.       */
  1133.       for (k=j; k<j2; ++k)
  1134.     partition[k]=i;
  1135.       numlines_l[i]=(j2-j);
  1136.       j = j2;
  1137.       if (j > BLKSIZE/2) break;
  1138.     }
  1139.   *npart_l_orig = i;
  1140.  
  1141.   /* compute which partition bands are in which scalefactor bands */
  1142.   { int i1,i2,sfb,start,end;
  1143.     FLOAT8 freq1,freq2;
  1144.     for ( sfb = 0; sfb < SBMAX_l; sfb++ ) {
  1145.       start = gfc->scalefac_band.l[ sfb ];
  1146.       end   = gfc->scalefac_band.l[ sfb+1 ];
  1147.       freq1 = sfreq*(start-.5)/(2*576);
  1148.       freq2 = sfreq*(end-1+.5)/(2*576);
  1149.              
  1150.       i1 = floor(.5 + BLKSIZE*freq1/sfreq);
  1151.       if (i1<0) i1=0;
  1152.       i2 = floor(.5 + BLKSIZE*freq2/sfreq);
  1153.       if (i2>BLKSIZE/2) i2=BLKSIZE/2;
  1154.  
  1155.       DEBUGF("old: (%i,%i)  new: (%i,%i) %i %i \n",bu_l[sfb],bo_l[sfb],
  1156.          partition[i1],partition[i2],i1,i2);
  1157.  
  1158.       w1_l[sfb]=.5;
  1159.       w2_l[sfb]=.5;
  1160.       bu_l[sfb]=partition[i1];
  1161.       bo_l[sfb]=partition[i2];
  1162.  
  1163.     }
  1164.   }
  1165.  
  1166. #endif
  1167.  
  1168.  
  1169.   /* compute bark value and ATH of each critical band */
  1170.   j=0;
  1171.   for(i=0;i<*npart_l_orig;i++)
  1172.     {
  1173.       FLOAT8 ji, freq, mval, bark;
  1174.       int k;
  1175.  
  1176.       ji = j;
  1177.       bark = freq2bark(sfreq*ji/BLKSIZE);
  1178.  
  1179.       ji = j + (numlines_l[i]-1);
  1180.       bark = .5*(bark + freq2bark(sfreq*ji/BLKSIZE));
  1181.       bval_l[i]=bark;
  1182.       //      j += numlines_l[i];
  1183.  
  1184.  
  1185.       gfc->ATH_partitionbands[i]=1e99;
  1186.       for (k=0; k < numlines_l[i]; ++k) {
  1187.     FLOAT8 freq = sfreq*j/(1000.0*BLKSIZE);
  1188.     assert( freq < 25 );
  1189.     //    freq = Min(.1,freq);    /* ignore ATH below 100hz */
  1190.     freq= ATHformula(freq);  
  1191.     freq += -20; /* scale to FFT units */
  1192.     freq = pow( 10.0, freq/10.0 );  /* convert from db -> energy */
  1193.     freq *= numlines_l[i];
  1194.     gfc->ATH_partitionbands[i]=Min(gfc->ATH_partitionbands[i],freq);
  1195.     ++j;
  1196.       }
  1197.  
  1198.  
  1199.       /*  minval formula needs work
  1200.       ji = j;
  1201.       freq = sfreq*ji/BLKSIZE;
  1202.       mval = Max(27.0 - freq/50.0, 0.0);
  1203.       //      mval = exp(-mval * LN_TO_LOG10);
  1204.  
  1205.       DEBUGF("%2i old minval=%f  new = %f  ",
  1206.          i,-log(minval[i])/LN_TO_LOG10,mval);
  1207.       if (mval < minval[i]) 
  1208.       DEBUGF("less masking than ISO tables \n");
  1209.       else
  1210.       DEBUGF("more masking than ISO tables \n");
  1211.       */
  1212.     }
  1213.  
  1214.  
  1215.  
  1216.   /************************************************************************
  1217.    * Now compute the spreading function, s[j][i], the value of the spread-*
  1218.    * ing function, centered at band j, for band i, store for later use    *
  1219.    ************************************************************************/
  1220.   /* i.e.: sum over j to spread into signal barkval=i  
  1221.      NOTE: i and j are used opposite as in the ISO docs */
  1222.   for(i=0;i<*npart_l_orig;i++)
  1223.     {
  1224.       FLOAT8 tempx,x,tempy,temp;
  1225.       for(j=0;j<*npart_l_orig;j++)
  1226.     {
  1227.           if (i>=j) tempx = (bval_l[i] - bval_l[j])*3.0;
  1228.       else    tempx = (bval_l[i] - bval_l[j])*1.5; 
  1229.  
  1230.       if(tempx>=0.5 && tempx<=2.5)
  1231.         {
  1232.           temp = tempx - 0.5;
  1233.           x = 8.0 * (temp*temp - 2.0 * temp);
  1234.         }
  1235.       else x = 0.0;
  1236.       tempx += 0.474;
  1237.       tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);
  1238.  
  1239. #ifdef NEWS3
  1240.       if (j>=i) tempy = (bval_l[j] - bval_l[i])*(-15);
  1241.       else    tempy = (bval_l[j] - bval_l[i])*25;
  1242.       x=0; 
  1243. #endif
  1244.       /*
  1245.       if ((i==cbmax/2)  && (fabs(bval_l[j] - bval_l[i])) < 3) {
  1246.         DEBUGF("bark=%f   x+tempy = %f  \n",bval_l[j] - bval_l[i],x+tempy);
  1247.       }
  1248.       */
  1249.  
  1250.       if (tempy <= -60.0) s3_l[i][j] = 0.0;
  1251.       else                s3_l[i][j] = exp( (x + tempy)*LN_TO_LOG10 ); 
  1252.     }
  1253.     }
  1254.  
  1255.  
  1256.   /************************************************************************/
  1257.   /* SHORT BLOCKS */
  1258.   /************************************************************************/
  1259.  
  1260. #ifdef NOTABLES
  1261.   /* compute numlines */
  1262.   j=0;
  1263.   for(i=0;i<CBANDS;i++)
  1264.     {
  1265.       FLOAT8 ji, bark1,bark2,delbark=.34;
  1266.       int k,j2;
  1267.  
  1268.       j2 = j;
  1269.       j2 = Min(j2,BLKSIZE_s/2);
  1270.       
  1271.       do {
  1272.     /* find smallest j2 >= j so that  (bark - bark_s[i-1]) < delbark */
  1273.     ji = j;
  1274.     bark1  = freq2bark(sfreq*ji/BLKSIZE_s);
  1275.     
  1276.     ++j2;
  1277.     ji = j2;
  1278.     bark2  = freq2bark(sfreq*ji/BLKSIZE_s);
  1279.  
  1280.       } while ((bark2 - bark1) < delbark  && j2<=BLKSIZE_s/2);
  1281.  
  1282.       /*
  1283.       DEBUGF("%i old n=%i  %f old numlines:  %i   new=%i (%i,%i) (%f,%f) \n",
  1284. i,*npart_s_orig,freq,numlines_s[i],j2-j,j,j2-1,bark1,bark2);
  1285.       */
  1286.       for (k=j; k<j2; ++k)
  1287.     partition[k]=i;
  1288.       numlines_s[i]=(j2-j);
  1289.       j = j2;
  1290.       if (j > BLKSIZE_s/2) break;
  1291.     }
  1292.   *npart_s_orig = i;
  1293.  
  1294.   /* compute which partition bands are in which scalefactor bands */
  1295.   { int i1,i2,sfb,start,end;
  1296.     FLOAT8 freq1,freq2;
  1297.     for ( sfb = 0; sfb < SBMAX_s; sfb++ ) {
  1298.       start = gfc->scalefac_band.s[ sfb ];
  1299.       end   = gfc->scalefac_band.s[ sfb+1 ];
  1300.       freq1 = sfreq*(start-.5)/(2*192);
  1301.       freq2 = sfreq*(end-1+.5)/(2*192);
  1302.              
  1303.       i1 = floor(.5 + BLKSIZE_s*freq1/sfreq);
  1304.       if (i1<0) i1=0;
  1305.       i2 = floor(.5 + BLKSIZE_s*freq2/sfreq);
  1306.       if (i2>BLKSIZE_s/2) i2=BLKSIZE_s/2;
  1307.  
  1308.       DEBUGF("old: (%i,%i)  new: (%i,%i) %i %i \n",bu_s[sfb],bo_s[sfb],
  1309.          partition[i1],partition[i2],i1,i2);
  1310.  
  1311.       w1_s[sfb]=.5;
  1312.       w2_s[sfb]=.5;
  1313.       bu_s[sfb]=partition[i1];
  1314.       bo_s[sfb]=partition[i2];
  1315.  
  1316.     }
  1317.   }
  1318.  
  1319. #endif
  1320.  
  1321.  
  1322.  
  1323.   /* compute bark values of each critical band */
  1324.   j = 0;
  1325.   for(i=0;i<*npart_s_orig;i++)
  1326.     {
  1327.       FLOAT8 ji, bark, freq, snr;
  1328.       ji = j;
  1329.       bark = freq2bark(sfreq*ji/BLKSIZE_s);
  1330.  
  1331.       ji = j + numlines_s[i] -1;
  1332.       bark = .5*(bark+freq2bark(sfreq*ji/BLKSIZE_s));
  1333.       /*
  1334.       DEBUGF("%i %i bval_s = %f  %f  numlines=%i  formula=%f \n",i,j,bval_s[i],freq,numlines_s[i],bark);
  1335.       */
  1336.       bval_s[i]=bark;
  1337.       j += numlines_s[i];
  1338.  
  1339.       /* SNR formula needs work 
  1340.       ji = j;
  1341.       freq = sfreq*ji/BLKSIZE;
  1342.       snr  = .14 + freq/80000;
  1343.       DEBUGF("%2i old SNR=%f  new = %f \n ",
  1344.          i,SNR[i],snr);
  1345.       */
  1346.     }
  1347.  
  1348.  
  1349.  
  1350.  
  1351.  
  1352.  
  1353.   /************************************************************************
  1354.    * Now compute the spreading function, s[j][i], the value of the spread-*
  1355.    * ing function, centered at band j, for band i, store for later use    *
  1356.    ************************************************************************/
  1357.   for(i=0;i<*npart_s_orig;i++)
  1358.     {
  1359.       FLOAT8 tempx,x,tempy,temp;
  1360.       for(j=0;j<*npart_s_orig;j++)
  1361.     {
  1362.           if (i>=j) tempx = (bval_s[i] - bval_s[j])*3.0;
  1363.       else    tempx = (bval_s[i] - bval_s[j])*1.5; 
  1364.  
  1365.       if(tempx>=0.5 && tempx<=2.5)
  1366.         {
  1367.           temp = tempx - 0.5;
  1368.           x = 8.0 * (temp*temp - 2.0 * temp);
  1369.         }
  1370.       else x = 0.0;
  1371.       tempx += 0.474;
  1372.       tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);
  1373. #ifdef NEWS3
  1374.       if (j>=i) tempy = (bval_s[j] - bval_s[i])*(-15);
  1375.       else    tempy = (bval_s[j] - bval_s[i])*25;
  1376.       x=0; 
  1377. #endif
  1378.       if (tempy <= -60.0) s3_s[i][j] = 0.0;
  1379.       else                s3_s[i][j] = exp( (x + tempy)*LN_TO_LOG10 );
  1380.     }
  1381.     }
  1382.  
  1383.  
  1384.  
  1385.  
  1386.   /* compute: */
  1387.   /* npart_l_orig   = number of partition bands before convolution */
  1388.   /* npart_l  = number of partition bands after convolution */
  1389.  
  1390.   assert(*npart_l_orig <=CBANDS);
  1391.   assert(*npart_s_orig<=CBANDS);
  1392.  
  1393.   
  1394.   *npart_l=bo_l[SBPSY_l-1]+1;
  1395.   *npart_s=bo_s[SBPSY_s-1]+1;
  1396.   
  1397.   /* if npart_l = npart_l_orig + 1, we can fix that below.  else: */
  1398.   assert(*npart_l <= *npart_l_orig+1);
  1399.   assert(*npart_s <= *npart_s_orig+1);
  1400.  
  1401.  
  1402.   /* MPEG2 tables are screwed up 
  1403.    * the mapping from paritition bands to scalefactor bands will use
  1404.    * more paritition bands than we have.  
  1405.    * So we will not compute these fictitious partition bands by reducing
  1406.    * npart_l below.  */
  1407.   if (*npart_l > *npart_l_orig) {
  1408.     *npart_l=*npart_l_orig;
  1409.     bo_l[SBPSY_l-1]=(*npart_l)-1;
  1410.     w2_l[SBPSY_l-1]=1.0;
  1411.   }
  1412.  
  1413.   if (*npart_s > *npart_s_orig) {
  1414.     *npart_s=*npart_s_orig;
  1415.     bo_s[SBPSY_s-1]=(*npart_s)-1;
  1416.     w2_s[SBPSY_s-1]=1.0;
  1417.   }
  1418.  
  1419.  
  1420.     /* setup stereo demasking thresholds */
  1421.     /* formula reverse enginerred from plot in paper */
  1422.     for ( i = 0; i < SBPSY_s; i++ ) {
  1423.       FLOAT8 arg,mld;
  1424.       arg = freq2bark(sfreq*gfc->scalefac_band.s[i]/(2*192));
  1425.       arg = (Min(arg, 15.5)/15.5);
  1426.  
  1427.       mld = 1.25*(1-cos(PI*arg))-2.5;
  1428.       gfc->mld_s[i] = pow(10.0,mld);
  1429.     }
  1430.     for ( i = 0; i < SBPSY_l; i++ ) {
  1431.       FLOAT8 arg,mld;
  1432.       arg = freq2bark(sfreq*gfc->scalefac_band.l[i]/(2*576));
  1433.       arg = (Min(arg, 15.5)/15.5);
  1434.  
  1435.       mld = 1.25*(1-cos(PI*arg))-2.5;
  1436.       gfc->mld_l[i] = pow(10.0,mld);
  1437.     }
  1438.  
  1439.   
  1440.   return 0;
  1441. }
  1442.